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Using an expansion in order parameters, the equation of motion for 
the centroid of globally coupled oscillators with natural frequencies taken 
from a distribution is obtained for the case of high coupling, low dispersion 
of natural frequencies and any number of oscillators. To the first order, 
the system can be approximated by a set of four equations, where the cen- 
troid is coupled with a second macroscopic variable, which describes the 
dynamics of the oscillators around their average. This gives rise to collec- 
tive effects that suggest experiments aimed at measuring the parameters 
of the population. 

Populations of coupled oscillators have been widely used as a framework 
for studying collective phenomena. Applications include chemical systems (e.g., 
reaction-diffusion systems in the oscillatory regime [1]), electronic circuits [2], 
and biological systems (e.g., spiking activity in neurons [3], synchronization in 
/3-cells [4], and glycolytic oscillations [5]). 

Much of the interest has been devoted to properties arising at the population 
level, especially in relation with synchronization phenomena [6]. Borrowing 
from statistical mechanics, such properties may be quantified in terms of order 
parameters [7], such as the mean value or centroid, for which phase diagrams 
can be constructed. 

In this Letter, we consider a population of N Hopf normal forms with a 
mean-field coupling: 



where the natural frequencies uij are taken from a given distribution and K 
denotes the coupling strength. 

Previous analyses [8] of system ([j]) (or of its phase reduction) have been 
restricted mainly to properties arising in the asymptotic regimes. In this case, 
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an insight into the qualitative behavior of the system has been obtained by 
observing the asymptotic properties of the centroid Z = (zj) (using the notation 

(fj) = jf fj)- 

Our aim is to describe the dynamics of system (|l|) by expanding it in order 
parameters or macroscopic variables (as they shall be called from now on in 
order to avoid confusion between control parameters and order parameters). 
The explicit deduction of their equations of motion shows how properties emerge 
from the microscopic level, and allows us to describe the dynamics of the system 
outside the asymptotic regime, a relevant issue in systems interacting with the 
environment or in experiments using perturbation response techniques. 

The macroscopic equations are derived through a change of coordinates from 
the microscopic variables Zj to macroscopic variables (like the centroid) that are 
averaged over the population. The main advantage of this procedure is that the 
macroscopic variables can be organized in a hierarchy, retaining, if the cou- 
pling is strong and the natural frequency distribution is narrow, only the lowest 
order terms of a series expansion. In this region of the parameter space, the 
equation for the centroid has a functional form reminiscent of Eq. (|l|), but it 
is also coupled with a second macroscopic variable that describes the dynamics 
of the oscillators around their average. As a consequence, with respect to the 
description of the system as a single macroscopic oscillator, a correction for the 
collective oscillations amplitude appears and a critical macroscopic perturba- 
tion (i.e., quenching, see later) can be identified. Analytical relations for such 
quantities are given and suggest simple experiments for determining, by means 
of purely macroscopic measures, the coupling strength and the variance of the 
natural frequencies distribution in real systems. Our results are valid for any 
number of oscillators, and for any (narrow) distribution of natural frequencies 
and are compared with those obtained by numerically integrating the original 
system (0). 

Deduction. Let us start by writing the positions of the oscillators in terms 
of their distance from the centroid: Zj = Z + ej. This expansion is useful 
when the coupling is strong, since in this case the dynamics leads the system to 
collapse on a configuration, the so-called locked state [9] , that is peaked around 
the centroid. In particular, the displacements from the centroid in the locked 
state converge to zero when the coupling strength is increased^. Moreover, when 
the displacements are wider than the value at the equilibrium, the dynamics is 
contracting, and thus the locked state furnishes an upper limit to the broadness 
of the initial configuration required for the approximations to hold. 

Eq. (|l]) now reads: 



fly 

-±= (1 - \Zf + iuj) Z + e 3 (iuj + 1-K- 2\Z\ 2 ) 

-Z*e*+o{\etf). (2) 

We then consider the time evolution of the centroid, differentiating its defi- 
nition: Z = (zj) and using Eq. (||): 

dZ/dt = d{zj)/dt = (dzj/dt) = (1 - \Z\ 2 + i(uj 3 ))Z + i(ujej) + (1 - K - 
2\Z\ 2 ){e ] )-Z\e*) + {o{\e J \ 2 )). 

By definition of a centroid, (e^) = (e*) = 0, and dZ/dt reduces to: 
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^ = (1 - \Z\ 2 + iuJ Q )Z + i( Uj ej) + o | 2 )) (3) 

where (w,*) = c^o is the average natural frequency. A zeroth-order expansion 
thus leads to the equation: 



— = (l-\Z\ 2 + t u; a )Z. (4) 

This expansion has the same functional form as the individual, uncoupled 
elements. It exactly describes the case of a population of oscillators with the 
same natural frequency ujq and with tj =0 Vj, the last condition being fulfilled 
if all the oscillators are assigned the same initial condition. In this (trivial) case, 
it shows the existence of a limit cycle of radius 1 and frequency ujq, and of an 
unstable focus in Z = 0. The equation is independent of the coupling strength, 
since the coincidence of the oscillators and the centroid is maintained in time. 
For the description of the parameter mismatch to be included, the term of lower 
order that must be then taken into account is: W :— (u)j ej). Having in this 
way defined a new macroscopic variable, we now derive its equation of motion. 

Since dW/dt = d{ujj€j)/dt — (u)jdzj/dt) — (ujj)dZ/dt, using the definition of 
W and Eq. (||), we obtain: 

^ = i ««;?) - u 2 ) Z + (1 - K - 2\Z\ 2 + iu ) W 

-Z 2 W* +i((cj j -u; ) 2 e j )+o({\e]\)). (5) 

The equation is closed with respect to Z and W, if again the higher order 
terms in the displacements are discarded, when the dispersion of the natural 
frequencies is sufficiently small. Since we are dealing with narrow distributions, 
the second degree term in the frequencies will be neglected. Note that this ap- 
proximation does not depend on specific symmetries or shapes for the frequency 
distribution. Calling a 2 — [uj 2 ) — uj 2 its variance, and combining Eq. (^|) with 
Eq. (ft), we get the following description at first order: 



«£ = {i-\z\ 2 + l u«)z + l w 

=ia 2 Z+ (l-K-2\Z\ 2 + iuj )W~Z 2 W*. 1 ' 

This first-order expansion shows explicitly how the centroid is coupled with 
the dynamics of the oscillators around it. This "internal" dynamics is quantified 
in the W variable, by means of which the parameter mismatch (through the 
variance of natural frequencies) and the coupling strength affect the system. 
These equations describe the dynamics of the centroid from initial conditions 
with small ej (and thus small W). Since in their derivation we have requested 
to have asymptotically a locked state, we will restrict our analysis to the case 
in which the system has an unstable focus surrounded by a limit cycle, which 
occurs if K > 1 + a 2 . 

Although we stop at this order, the deduction can be carried further, and 
macroscopic descriptions at higher orders derived, applying the same method: 
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the terms that have been discarded in Eqs. (||) and (||) define new macroscopic 
variables, and their equations of motion can be derived by differentiating their 
definition. 

Analysis. If T^T > 1 + a 2 and ej are small, it is easy to see that the behavior 
of the system is essentially described by the dynamics on the invariant and at- 
tracting manifold on which the two macroscopic variables are orthogonal. This 
condition is fulfilled, and the description exact, when the natural frequency dis- 
tribution and initial configuration are both symmetric. Setting Z = R e i (0+'"o*) 
and W = we l ( e+uot \ Eqs. (||) reduce on = 6+ir/2 to the amplitude equations: 



--(l-R 2 )R~w 
-- a 2 R + (1 - K - R 2 ) w. 

System (Q) can be analyzed with planar methods and furnishes both an 
estimate for the amplitude of the locked oscillations and a description of their 
transient behavior. 

Let us first find the equilibrium values for R = \Z\ and w = \W\. Setting 
a = K/2 — y/ (K/2) 2 — a 2 and using Eq. it is easy to show that W and 
Z display synchronous oscillations (the second with a phase delay of n/2) of 
amplitudes: 



Ri = vl — a, w\ — ol\/\ — a. (8) 

Taking into account that in our approximation the ratio a 2 /K is small, this 
expression approximately results in: 

When these expressions are compared with the zeroth-order amplitudes Rq = 
1 and wo = 0, it appears that the term a 2 /2K gives the first order correction to 
the amplitude of the oscillations due to the mismatch of the natural frequencies. 
It is worth remarking that previously derived formulae for the amplitude of Z 
were only implicitly related to the parameters of the population, through a 
self-consistency integral [9], and in the limit of large N. 

Let us now address the structure of the phase space. This is particularly 
relevant when the perturbation response is addressed, since it qualitatively and 
quantitatively determines the features of the transient after the system is ini- 
tialized in a configuration out of equilibrium. A fundamental example is the 
quenching technique [10], which consists in damping the oscillations of the sys- 
tem by displacing it onto the stable manifold of the unstable focus. This can 
be experimentally observed as a long-term vanishment of the oscillations am- 
plitude. While in the case of the zeroth-order expansion the stable manifold 
reduces to the focus itself, in the case of the first-order expansion it consists of 
a two-dimensional variety, which on the invariant plane reduces to the stable 
eigenspace of the unstable focus: 



w-(K-a)R= 0. (10) 
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In quenching experiments the system lies in the locked state and is then 
macroscopically perturbed through a rigid displacement of the whole popula- 
tion. Since only the amplitude of the collective oscillations is measured, the 
critical displacement is identified by means of the quenching radius, that indi- 
cates the amplitude of the initial state for which the oscillations are suppressed. 
According to Eq. (0), this radius is given by: 



_ aVl-a _ a 2 

Again, we notice that the expression reduces to the zeroth-order value, that 
is zero, in the limit of vanishing a and large K . 

Moreover, as illustrated in Fig. [l], the structure of the phase space allows 
us to predict three qualitatively different kinds of transient amplitude increase. 
In one case (I), the system behaves like an ordinary oscillator, the impulsive 
perturbation being followed by a monotonic increase in amplitude. In the second 
case (II), the amplitude first decreases, while the centroid is approaching the 
origin, and then increases again, leading at the same asymptotic solution as in 
the previous case. Finally (III), if the initial configuration width is large enough, 
it will first reduce to zero and then increase again, but the asymptotic solution 
will be in phase opposition with respect to the initial oscillations. 

Comparison. Let us now compare the behavior of the equations for the 
macroscopic variables Eq. (^) with the numerical simulations of the original 
system described by Eq. (Q). Fig. ^ compares the transient behavior of the full 
system (initialized in three configurations having the same average Z , but dif- 
ferent W) with the one predicted according to Eq. (||) . The transient behaviors 
shown in Fig. [I] are actually observed: the accuracy of the approximation holds 
along the whole trajectory and remains nearly unchanged for any population 
having the same a and K. Fig. | compares the estimated and numerically 
computed values of the asymptotic amplitude of the collective oscillations and 
the quenching radius. 

We considered a population of globally and strongly coupled oscillators with 
narrow natural frequency distribution and showed how the dynamics can be re- 
duced, via an expansion in order parameters, to a low-order system containing 
all the essential information (for the first-order expansion, two complex vari- 
ables instead of N and two parameters instead of N + 1). Qualitative as well as 
quantitative results were then given for collective properties of experimental rel- 
evance. Although we restricted our attention to simple oscillators, the presented 
method can in principle be applied to other populations of globally and strongly 
coupled elements with small parameter mismatch. Moreover, the macroscopic 
description may furnish an instrument for studying collective behaviors outside 
the phase locking region. In particular, the number of the relevant terms in 
the expansion is expected to change as the system approaches the bifurcation 
boundaries, linking the onset of new regimes to a change in the dimensionality 
of the macroscopic system. 
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Figure 1: Phase portrait of the reduced system Eq. (^|). The equilibria (E\ and 
E2) and the eigenspaces of the origin are indicated. The system is initialized in 
three points (circles) giving rise to qualitatively different trajectories. 
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Figure 2: The transient behavior predicted by Eq. @ (solid line) is compared to 
that of the full system Eq. ([[]) (triangles) and of its zeroth-order approximation 
Eq. ([|) (dotted line) for a 2 — 0.5 and K = 3. Three initial states are chosen, 
having the same centroid's position \Z\, but different configuration, and thus 
different \W\. Populations with different size and frequency distribution are 
considered: N = 800, Gaussian distribution (A); N = 800 uniform distribution 
(V); N = 5, uniform distribution (<); N = 2 (t>). 
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Figure 3: The estimated values for the amplitude of the centroid's oscillations 
(Eq. |^) and the quenching radius (Eq. [Tl]) versus the coupling constant K (solid 
lines) are compared to the numerically computed ones (triangles), for a 2 = 0.2 
(cases a and b) and a 2 — 0.5 (cases c and d). The accordance is higher for more 
concentrated distributions, but is again almost unchanged when populations of 
different size and frequency distribution are considered (symbols as in Fig. ^|). 
The accuracy increases with the coupling, approaching at the meantime the 
zeroth-order values \Z\ = 1 and R q — 0. 
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